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We incorporate a massless scalar field into a 3-dimensional code for the characteristic evolution 
of the gravitational field. The extended 3-dimensional code for the Einstein-Klein-Gordon system 
is calibrated to be second order convergent. It provides an accurate calculation of the gravitational 
and scalar radiation at infinity. As an application, we simulate the fully nonlinear evolution of an 
asymmetric scalar pulse of ingoing radiation propagating toward an interior Schwarzschild black 
hole and compute the backscattered scalar and gravitational outgoing radiation patterns. The 
amplitudes of the scalar and gravitational outgoing radiation modes exhibit the predicted power 
law scaling with respect to the amplitude of the initial data. For the scattering of an axisymmetric 
scalar field, the final ring down matches the complex frequency calculated perturbatively for the 
1 = 2 quasinormal mode. 

PACS numbers: 04.40.-b, 04.25. Dm, 04.40.Nr, 04.40. Dg 



Numerical relativity has made significant advances in simulations which will be helpful toward achieving the final 
goal of providing useful information for the detection of gravitational waves, although the goal of achieving long term 
simulations of binary black holes has not yet been met. At some level, we have entered an era in which Einstein's 
equations can effectively be considered solved at the local level; i.e. there are several codes which provide accurate 
evolutions in a sufficiently small space-time region where the gravitational field has small variation when expressed in 
a Riemann normal coordinate system. However, at the global level, the computation of gravitational radiation from 
black hole space-times remains a highly challenging problem, in which there are vastly different time scales at play in 
the inner and outer regions and in which the proper choice of gauge is far from obvious. 

Most work in numerical relativity is based upon the Cauchy 3 -f 1 formalism (see, for instance 1]). A different 
approach, which is specifically tailored to study radiation, is based upon the characteristic initial value problem. This 
approach, which has been successful in computing gravitational radiation in single black hole space-times, is the focus 
of attention in this paper. 

Bondi, Sachs and Penrose have proposed formalisms to deal with radiation which today are the cornerstones 

for the characteristic numerical formulation of general relativity. All schemes for characteristic evolution have a 
common structure. The main ingredient is the space-time foliation by null hypersurfaces u = constant, generated by 
a set of bi-characteristic null rays , with a coordinate A varying along them. In null coordinates [u, A, x^) the 
main set of Einstein's equations has the form 



where F represents variables intrinsic to a single null hypersurface u, G represents the evolution variables and. Hp 
and Hq are nonlinear operators intrinsic to the null hypersurface u. Besides these main equations, there is a set of 
constraint equations which are satisfied via the Bianchi identities if they are satisfied on an inner boundary. At null 
infinity, these constraints can be interpreted as conservation conditions governing energy and angular momentum. 

The numerical implementation of the characteristic method consists in evolving a given field (e.g. scalar, electro- 
magnetic or gravitational) on a family of null hypersurfaces along a discrete sequence of retarded time steps. A stable, 
second order accurate, fully nonlinear, 3-dimensional code (the PITT code) has been based upon this characteristic 
formalism. Its implementation, tests and results have been presented in a series of papers [E S S • The code 
poses data on an initial null hypersurface and on an inner worldtube boundary, and evolves the exterior space-time 
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out to a compactified version of null infinity. The code calculates waveforms at null infinity 0,^3; tracks a dynamical 
black hole and excises its internal singularity from the computational grid 0, 0| . 

The PITT code uses an explicit finite diiTerence evolution algorithm based upon retarded time steps on a uniform 
3-dimensional null coordinate grid |^ ^3 ■ crucial ingredient of the code is the S~module which incorporates a 
computational version of the Newman-Penrose 9- formalism 13] . The 9-module covers the sphere with two overlapping 
stereographic coordinate grids (North and South). It provides everywhere regular, second order accurate, finite 
difference expressions for tensor fields on the sphere and their covariant derivatives. The 9-calculus simplifies the 
underlying equations, avoids spurious coordinate singularities and allows accurate differentiation of tensor fields on 
the sphere in a computationally efficient and clean way. The approach leads to a straightforward implementation of 
a second order convergent, finite difference code. The PITT code has undergone recent improvements to increase the 
accuracy in the simulation of highly distorted single black hole space-times |p| and in the calculation of the Bondi 
news function which describes the radiated waveform at null infinity X"*" [T3. Il5j|. The characteristic code has been 
extended to include a naive hydrodynamical treatment, which has been applied to model a fluid ball falling into a 
Schwarzschild black hole and to model a polytrope in orbit near a Schwarzschild black hole 0| . 

Our work here is motivated by the possibility of using characteristic evolution to model bosonic stars in orbit 
about a black hole. Besides the possible importance of such stars astrophysical and cosmological aspects [T8llT9ll20l |. 
their existence as solitonic solutions [2lL \22L |23| offers a way to model a compact object without the difficulties of a 
hydrodynamic description of the matter field. The bosonic field can be evolved using the same characteristic techniques 
developed for the gravitational field. With this goal in mind, here we incorporate a scalar field in the PITT code. We 
demonstrate that the code provides a stable, convergent evolution of a 3-dimcnsional massless Einstein-Klein-Gordon 
(EKG) system in the region exterior to a black hole. 

This extension also allows the study of the nonlinear interaction between gravitational and scalar radiation in a 
fully 3-dimensional black hole space-time. Previous studies have measured the nonlinear effects in the scattering 
of a gravitational wave off a black hole, using codes based upon characteristic evolution |^, and upon Cauchy 
evolution In the purely gravitational case, the spherical harmonic modes of the gravitational waves are 

mixed by quadratic (and higher order) interactions. In the case of the massless EKG system, the scalar wave modes 
are only mixed indirectly through their quadratic coupling (via their stress energy tensor) to the gravitational field. 
As a result, scalar wave modes excite gravitational wave modes at the quadratic order but they excite other scalar 
modes only at a cubic order. Thus the scalar-scalar interaction is weaker and less interesting physically than the 
scalar-gravitational interaction. The production of gravitational waves by the scattering of scalar waves has many 
mathematical features in common with the production of gravitational waves by the motion of fluid bodies. The 
results presented here demonstrate the versatility of the characteristic code to compute gravitational waves generated 
by scalar sources and are an encouraging step toward eventually treating a bosonic star in orbit about a black hole. 

The massless scalar field coupled minimally with gravitation has been thoroughly studied in ID. Choptuik '2^12^ 
carried out the first numerical study of critical behavior in the collapse of massless scalar fields in the context of 
spherical symmetry. Brady, Chambers and Goncalyes [2^ studied critical behavior in the collapse of spherically 
symmetric massive scalar fields. Seidel and Suen pOl studied the formation of bosonic stars from massive, complex 
scalar fields. Balakrishna, Seidel and Suen js^ evolved self-gravitating massive scalar fields to study the stability 
of bosonic star configurations. Recently, axisymmetric simulations of of a complex scalar field have been used to 
study angular momentum in critical phenomena 'si'l and fully three dimensional simulations have been used to study 
bosonic stars (32] and gravitational collapse ,33j . 

A code based upon both incoming and outgoing null cones has been used in a combined global treatment of 
future infinity and a black hole horizon for an Einstein- Klein-Gordon field with spherical symmetry 0. Marsa and 
Choptuik [3j| used Eddington-Finkelstein coordinates to study the dispersion of scalar waves in ID. Siebel, Font and 
Papadopoulos j35l ] studied the interaction between a massless scalar field and a neutron star modeled as a perfect 
fluid. The first characteristic code in Bondi coordinates for the self-gravitatin g sc alar wave problem was constructed 
by Gomez and Winicour [s^. Subsequently Gomez, Schmidt and Winicour js^] applied the characteristic code to 
study the radiation tail decay of a scalar field. Barreto et al. '3^ also used this characteristic code to study the 
instability of a topological kink in the configuration of the scalar field. Recently, 2-dimensional and 3-dimensional 
Cauchy simulations of the massless scalar field have been carried out in the perturbative and linearized regimes |39l l40| ] 
(and references therein). 

The paper is organized as follows. In Sec. ^ we give the field equations in the characteristic form for a massive, self- 
interacting, complex scalar field minimally coupled with gravity. In Sec. lIIII we describe the numerical implementation. 
In Sec. II VI we present convergence and stability tests and results for the scattering of a scalar field off a black hole. 
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II. FIELD EQUATIONS 

We use coordinates based upon a family of outgoing null hypcrsurfaces. We let u label these hypersurfaces; 
{A = 2,3) label the null rays and r be a surface area coordinate. In the resulting x" — {u,r,x^) coordinates, the 
metric takes the Bondi-Sachs form (3, Q 

ds^ = -[e^'^il + W/r) - r^hABU^U^]du^ - 2e^'^dudr 

-2r^hABU^dudx'^ + r^hAsdx'^dx^ , (3) 

where W is related to the more usual Bondi-Sachs variable V hy V = r -\- W, and where h^^hsc — and 
det{hAB) — det{qAB), with qab a unit sphere metric. We also use the intermediate variable 

Qa = r^e-^f'hABU^- (4) 

We work in stereographic coordinates x"^ = {q,p) in which the unit sphere metric is 

QABdx^dx^ = -^{dq^ + dp^), (5) 

where 

P=l+p^ + q^. (6) 

We also introduce a complex dyad qA defined by 

q^ = ^ih^), (7) 

with i = ^—1. For an arbitrary Bondi-Sachs metric, Hab can then be represented by its dyad component 

J = /iAs9^g^/2, (8) 

with the spherically symmetric case given by J = 0. The full nonlinear Hab is uniquely determined by J, since the 
determinant condition implies that the remaining dyad component 

K = hABq'^f/2 (9) 

satisfies 1 — — J J. We also introduce spin-weighted fields 

u = u\a, Q = Q^g^, (10) 

as well as the (complex differential) operators 3 and 9 (see jQ] for full details). 

The null cone problem is normally formulated in the region of space-time between a timelike or null world tube T 
and X+, with initial data J given on an initial null cone u = 0. Boundary data for the metric variables /3, Q, U , W 
and J arc required on T. We represent on a finite grid by using a compactified radial coordinate x — t'/{1 + r), 
in terms of which all code variables are globally regular. 

The general field equations for a complex scalar field coupled minimally to gravity are 

Rab = &Tr{Tab - ^9abT), (11) 



with 



Tab = ■^{4>,a4>,b + 4>,a4>,b) - gab {^9"'^ 4> ^dp 4 + (12) 



where -F(0, (p) includes a possible mass term and self-interacting potential. The scalar field obeys the wave equation 

dF 

□0-^=0. (13) 

In this paper we treat the case of a massless, real scalar field where F = 0. The field equations then reduce to 

Rab = Sn(j)^a(t),b- (14) 
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The corresponding Bondi-Sachs hypersurface equations are 

= 2^r(0,,)2 + Np, (15) 

{r'^Q)^r = l&T^r^(j},r^(j} - r2(§J + ^K)^r + 2r^^{^)^r + Nq, (16) 

U,r^ —{KQ~ JQ), (17) 

le^^n - 1 - e'^SBe'^ + -^\r*(dU + W)] r + Nw, (18) 
2 



where 



n = 2K~ BBK + i(g2 J + g2 J) + -^(gjgj - gjgj); (19) 

2 4iv 



and the evolution equation for J is 



2(rJ),„, - [F/r(rJ),,]^, 



— e2^(g0)2 - -{r^dU) r + -e^B^e" - (W/r) , J + Nj. (20) 
r r r 

The expressions for A'^, Nq, Nw and A^j are given in 10]. The wave equation 1)1 3|1 for the scalar field is 

2(r0),„, - [V/r{rcj,)^rir = -{W/r)^r^ + N^, (21) 



where 



e2/3 



iV^ = — (iV^i - N4.2 + N^s) - 2^04 - N^5, (22) 



AT^i = K{dd(p + B(3d(l> + dpdcj)) , (23) 



N^2 = ^[3J3</' + SJ^<I) + 2{JdpB<j) + Jdf3B(j)) + JB^cp + jg^,?^], (24) 



= + jgjg0 + Jgjg0 + jgjg^), (25) 



N^i = (j)^r{W + W) + 2{m(t>,r + UQ(l),r) + U^rQ4> + U,rQ(t), (26) 



iV^5 = Ud(b + Udcb. (27) 

The data required on the initial null cone are J and 0, which constitute the evolution variables. The remaining 
auxiliary variables can then be determined on the initial null cone by explicit radial integration in the following order: 
f3 from Eq. (fT^ . Q from Eq. (fT?)|l . U from Eq. lfT7)l and W from Eq. ifT^ . The evolution equations (jSOJ and if^ can 
then be used to find J and </> (in that order) on the "next" null cone. Boundary data on F is required for the scalar 
field (j) in addition to the gravitational variables J, /3, Q, U, and W. 

In concluding this section we point out that we obtain the Bondi news function by carrying out a transformation 
to an inertial frame at X+ . The news is then expressed in terms of the two standard polarization modes A^_|_ and , 
as expressed in the computational g formalism (for details see Ref. 10]). 
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III. NUMERICAL IMPLEMENTATION 

To solve Eqs. (|15|) - (|21|) we follow the strategy developed for the vacuum and fluid matter cases 0,0] based upon 
a second-order accurate finite difference approximation. The details follow. 

A. Numerical grid 

We define a numerical grid with coordinates (u„, Xi, qj,Pk) — {nAu, l/2 + {i~l)Ax, — f + {j + 3)Aq, — f + {k + 3)Ap), 
where the spatial indexes range from i ~ 1 . . . N^, {j, k) = 1 . . . N(^, with 2Ax = l/{Nx — l), and Aq = Ap = 2/ {N(^ — 5). 
Using finite differences to discretize the equations, we center the derivatives at {n + 1/2, i — f/2, j, k). The evolution 
proceeds with time step Au, subject to a Courant-Friedrichs-Levy (CFL) condition 

B. Hypersurface equations 

The hypersurface equations (|f 5|I - Hf 8(1 are discretized as in with the right hand sides now including the scalar 
field evaluated at the mid-points Xj^_i on the radial grid, which is straightforward to interpolate from the scalar field 
values at the integral grid points Xi. 

C. Evolution equations 

The evolution equation (|20|) for J is treated exactly as in after modifying the right hand side to include the 
scalar field terms. Thus, we use a Crank-Nicholson scheme to solve the finite difference version of Eq. H2(J|) written in 
terms of the compactified radial coordinate x. This scheme introduces some numerical dissipation that stabilizes the 
code even in the regime of very high amplitude fields jlC|. 

We integrate the evolution equation for the scalar field using the null parallelogram marching algorithm [111 . 
Equation H21|) is recast in terms of the 2-dimensional wave operator 



corresponding to the line element 



= e-2/3[2(r0),,„ - ir-'Vir<t>),r).r] (28) 



da^ = 2l^^n^^dx''dx'' = e'^^[r~'^Vdu + 2dr], (29) 



where — u^f^ is the normal to the outgoing null cones and is an inward normal null vector to the spheres of 
constant r. The evolution equation for the scalar field then reduces to 

e"'3n(2)(r0) =7^, (30) 

where 

n = -{W/r),r(l) + N4,. (31) 

Because all 2-dimensional wave operators are conformally flat, with conformal- weight —2, we can apply to (|30|) a 
flat-space identity relating the values of rcf) at the corners P, Q, R and 5* of a null parallelogram A, with sides formed 
by incoming and outgoing radial characteristics. In terms of rcj), this relation leads to an integral form of the evolution 
equation for the scalar field: 



{rcp)Q = {rcp)p + {t4>)s - {rcp)R + dudrU. (32) 



The corners of the null parallelogram cannot be chosen to lie exactly on the grid because the velocity of light in terms 
of X coordinate is not constant. The values of r<j) at the vertices of the parallelogram are approximated to second-order 
accuracy by linear interpolations between nearest neighbor-grid points on the same outgoing characteristic. Then, 
by approximating the integrand by its value at the center C of the parallelogram, we have 

(T<t>)Q = {'r<t>)p + {r4))s - {r(j))R + ^Au{rQ -rp + rs - rpjUc- (33) 
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In order to apply this scheme globally we must also take into account technical problems concerning the order of 
accuracy for points near I"*". For this purpose, it is convenient to renormalize 1)32(1 by introducing the intermediate 
variable $ = {r(j)){l — x) — x4>. We choose (j> = initially at X+. The finite difference version of the evolution equation 
for the scalar field preserves this property. After this substitution, the evolution equation for <I> becomes 

$Q = ^xqAuTCc + f $p + ^xpAunc] + f $s + IxsAuHc] - f^n + ^xrAuHc] ■ (34) 

4^ 1-xp \ 4 J l-xs \ 4 J l-XR \ 4 / 



Since is defined by linear interpolation, we find (for 2 < i < Nx) 



and for z = 2 or z = Nr^ 



<i>QAx-^1+,\x,-XQ) 

Xi{xQ - Xi-^l) 



(35) 



(36) 



IV. SCATTERING OF A MASSLESS SCALAR FIELD OFF A SCHWARZSCHILD BLACK HOLE 

A. Setting the initial and boundary data 

We use the code to evolve an initial configuration where the gravitational data corresponds to a Schwarzschild black 
hole of mass M = 1 superimposed with initial data of compact radial support for the massless scalar field. The initial 
gravitational data are 

J{u = 0,x,x^) = 0. (37) 

In the absence of the scalar field, these data would generate the Schwarzschild solution throughout the exterior 
Kruskal quadrant. In the presence of the scalar field, the inner Schwarzschild white hole horizon is undisturbed but 
the black hole horizon is dynamically deformed. Note that there are no elliptic constraints in prescribing characteristic 
gravitational data, as opposed to the corresponding Cauchy problem. 
The data at the inner boundary F, located at r = 2M, are 

/3 = 0, Q = 0, C/ = 0, W^~2M, J = 0, = 0. (38) 

We take as the initial data for the scalar field 

,f n ^ ( Xir - ra)*{r - rb)*G{q,p) iire[Ra,Rb], /gn^ 

r0(u = O,r,g,p) = |Q otherwise. (^^^ 

Such data, prescribed here on an outgoing null hypersurface, corresponds to an ingoing pulse of scalar radiation. For 
G{q,p) — 1, the pulse would be a spherically symmetric shell with amplitude A located between ra and rt- For the 
convergence test in this paper, in both hemispheres we take G to be the function 

^('^'^^) = |o otherwise, ^^0) 

where = {q — qs)"^ + (p — Ps)^. Axisymmetric and non-axisymmetric initial massless scalar field configurations can 
be obtained depending on the values of qs^ Ps and fi. For q^ = Ps = 0, the initial profile of r(j) is axisymmetric about 
the poles and goes to zero at q^ + p^ = jjL. For other values of {qs,Ps), ^ is in general asymmetric. 

B. Monitoring convergence 

We measure convergence of the scalar field in terms of the norm 

2 



(41) 
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where the integration volume is taken over the entire null hypersurface exterior to r = 2M . We use the four-point 
formula 



f>'j±lM 



^J±l,fc±l)^ 



dq dp 



(42) 



for the contribution from cells contained completely in the upper and lower hemisphere. For a cell that crosses the 
equator, we approximate the equator by a straight line and evaluate the cell's additional contribution to H42|) by 



Q 



P 



P 



P 



],k±l 



P 



]±l,k±l 



dS. 



(43) 



where dSj^k is the trapezoidal surface element. The integration in x is straightforward; we use a second-order accurate 
Simpson's formula. The computation of l|41(l is the same order of finite difference approximation as used to solve 
(Cai-IIHl). 

For the convergence test we take Va = 3.5 and = 12, A = 10^^, /i = 0.3 and qg = Ps = 0. This gives an 
axisymmetric, equatorial reflection symmetric configuration corresponding to two ingoing pulses centered about the 
axis of symmetry. The following grids were used: 



• Coarse, Ux =41, n, 

• Medium, n. 



q 

81, na 



n.r. 



25 

= 45 



• Fine, Ux = 161, Uq ~ Up — 85, 

for which Ax and Aq — Ap scale as 4 : 2 : 1. Assuming that the quantity Q behaves as Q = a + 6A", it can be shown 
that the convergence rate is 



log2 



(44) 



where Qc, Qm and Qf refer to computed values of Q using the coarse, medium and fine grids, respectively [13 ■ The 
results in Table I show that the 3-dimensional EKG code is second order convergent in amplitude. 

We also measure the convergence in phase. It can be easily shown that the order of convergence in phase is expressed 
by 



log 



Qc 



where 



and 



Qmf 



Qmf 
{4>c - 4'm) 



(45) 



(0m - (t)f) 
P 



dqdpdx 



dqdpdx 



are calculated at the same grid points and at the same time by sub-sampling from the fine and medium grids to the 
coarse one. The results in Table II confirm that our code is also second order convergent in phase. 



C. Verifying the stability of the EKG-PITT code 



We have tested the code by carrying out numerical experiments which confirm that it is stable, subject to the 
CFL condition. The stability test consists of evolving random, initially localized scalar data of small amplitude for a 
sufficiently long time. The data is chosen to be random so that all possible numerical frequencies at a given resolution 
are present. Hence, any unstable mode is most likely to be excited and to dominate the evolution in the number of 
time steps taken. In particular, data of amplitude ~ 10~^ was run from m = to m = 50, corresponding to 2 x 10"* 
time steps for the coarse grid and 4 x 10"* time steps for the medium grid. During this time, the norm Q{u) initially 
decreases due to scattering of the scalar field. Figure ^ shows that the norm remains bounded at the end of the run, 
with a value of Q ~ 6 x 10"^''. 
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D. Scattering, scaling and mode coupling of axisymmetric pulses off a black hole 

We first consider initial data H39|l describing an £ = 2, m = massless scalar pulse propagating inward toward a 
Schwarzschild black hole, with the boundary conditions H38II and initial gravitational data H37|l . The parameters of the 
initial data are = 3, rf, = 5, A = 10~^. The radial location of the initial pulse is chosen for computational economy 
since nonlinear effects are weak for r >> 3M. Figure |2 shows the computed time dependence of rcj) at X"*", overlaid 
with the analytic curve for the corresponding quasi-normal mode calculated from perturbation theory lis","!^. After 
u = lOM, the ringdown of the computed waveform agrees with the complex frequency uj = 0.4875 — i0.098, obtained 
in the perturbative treatment of a scalar field on a Schwarzschild background. 

Next we consider the scattering of an ^ = 1, m = scalar pulse, with = 3, rt, = 5, for a range of amplitudes A. 
Figure 01 show the rescaled scalar waveform r0/A scattered to in the q ^ p = direction; and Fig. 01 shows the 
rescaled component /A^ of the gravitational news function at 1+ in the q — p — 0.5 direction. All the curves overlap 
for A = 10~^, lO""*, 10"'^. Thus it is manifest that rcj) has a linear dependence on amplitude and has a quadratic 
dependence, as expected. The results also show that £ = 1 is the lowest scalar mode that generates gravitational 
radiation. The quasinormal frequency calculated using perturbation theory is <I> = 0.292935 — lO. 09766 A 
comparison of the numerical and perturbative scalar waveforms over the time interval m < 20 reveals good agreement 
between the oscillation periods. However, during this time interval the numerical waveform does not yet display a 
clean quasinormal decay. Further investigation of this feature would require longer runs with higher resolution. 

Now we consider scattering of an axisymmetric scalar "blob" (|39|l of compact support < 0.4 on the north 

patch, with ra = 3, rb = 5, A = 10^^ and 

- 2 _ Q 4-|4 

The initial angular structure is dominated by the £ = 1, rn = harmonic. Figures [31 and (HI display snapshots of 
the evolution of rcj) and the spin- weight invariant J J at respectively, on the North hemisphere, obtained with a 
resolution of 89 x 89 angular grid points and 101 radial grid points. The radiation pattern has sharp angular features 
arising from the phase differences between the backscattering at different angles. The angular resolution is sufficient 
to reveal the main features of the evolution in fairly smooth snapshots. The numerical evolution of the axisymmetric 
initial data preserves the axisymmetry. The gravitational invariant J J vanishes at the pole g = p = 0, as required by 
axisymmetry. Figure [3 shows a global view of how the scalar field is backscattered to infinity while approaching the 
black hole. 

FigureslHlandOdisplay the rescaled radiation amplitudes r0/A and 7V+/A^ at 1+ produced on the North hemisphere 
by the scattering of the blob with varying amplitude A. Once again, the scalar field radiation pattern scales linearly 
and the gravitational polarization mode scales quadratically for incident amplitudes up to A = 10^^. Note in 
Fig. O that the scaling of 7V_|_ breaks down at the very small A = 10~^^ amplitude where the quadratic response is 
below numerical error. 

It is also possible to analyze the nonlinear content of the scalar waveform at 2^ by decomposing it in the spin-0 
spherical harmonic amplitudes 



Fi,n = <j> r4>Yi„,dVt. (47) 

These are computed by second order accurate integration over the sphere (with solid angle U, — 47r). We vary the 
amplitude of the ingoing axisymmetric blob from A = 10~^^ to A = 10~^ and carry out the simulations with a 
resolution of 25 x 25 angular grid points and 41 radial grid points (the coarse grid in the convergence tests). This 
is sufficient to illustrate the qualitative behavior with reasonable computational expense. We display the results for 
Fira{u)/\ with ? = 2, 4, 6 and TO = 0, in Figs. [1010 

The graph of F2o{u)/\ in Fig. llUI shows that the radiation amplitude varies linearly with input amplitude A, even 
at the larger amplitudes. The graphs in Fig's [TTl and 1121 of F4o(u)/A and i^6o(w)/A show similar behavior, i.e. a 
predominantly linear response, without nonlinear effects on the scalar field radiation. Thus it appears that up to an 
amplitude of A = 10~^ that the scattered waveform of the scalar field is in excellent agreement with results obtained 
by perturbation theory on a Schwarzschild background. This is in accord with expectations since the nonlinear effects 
are introduced by the gravitational field and thus are O(A^) compared with the 0(A) linear terms. 



E. Scattering of a non-axisymmetric pulse off a black hole 

In order to confirm the robustness of the implementation, we consider the scattering of an asymmetric scalar blob 
of compact support by a Schwarzschild black hole. The initial pulse (|39|) has parameters A — 10^^, = 3.5, r^ = 12, 
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/i = 0.03, with offset Qs = 0.2, ps = 0.3 (which brealcs the axisymmetry) . The evolution is carried out on a 85 x 85 x 101 
grid. 

Figures IT^ and ITU show surface plots at X+, at four different times, of the resulting scalar radiation rcf) and the 
gravitational wave energy flux |7Vp — \N^\'^ + lA^xP, respectively. As the evolution proceeds, the nonlinear effects 
produce a rich structure in the gravitational radiation. It is also clear from the snapshots of the scalar field that much 
finer angular grids are necessary to obtain smooth results in the asymmetric case. 

F. Speed of the code 

The tests were performed on Linux machines with a single AMD Athlon T-Bird processor running at 900MHz. 
One unit of physical time on a grid of 85 x 85 x 101 points takes about 21 hours. 

V. CONCLUSION 

In this paper we have incorporated a massless scalar field into a 3-dimensional characteristic code for the Einstein 
field equations. The code has been verified to be stable and convergent. A scalar field provides a clean way to simulate 
a variety of physically interesting scenarios: 

• Highly-perturbed black holes. The single black hole formed just after the merger of a binary black hole, or just 
after the asymmetric collapse of a star, would be highly distorted. A fully nonlinear general relativistic code is 
necessary to probe this regime. This requires a stable code that can deal with a single, dynamical, black hole 
and a "perturbing" agent, as can be modeled by a scalar field. Recent studies of gravitational waves incident on 
a black hole j24] have revealed interesting phenomena arising from such highly nonlinear perturbations of the 
gravitational field. The scalar field provides an alternative perturbation to shed light on the robustness of these 
effects. 

• Toroidal distributions of "sources" around black holes. Toroidal distributions of matter fields around black holes 
play a major role in models describing active galactic nuclei, quasars and even gravitational wave sources. A 
successful simulation of these systems will require, at the very least, robust general relativistic codes appropri- 
ately coupled to relativistic hydrodynamics. A toroidal scalar field distribution around a black hole represents 
a intermediate step for assessing the performance of the implementation. Although this has limitations, it does 
provide a geometrical set up similar to the one expected in real astrophysical sources. Therefore, it can shed 
light on issues such as the gravitational wave output of physical systems, the characteristics of the expected 
waveforms, the radiated angular momentum, etc. Additionally, the equivalence of a massless scalar field to an ir- 
rotational stiff fiuid provides a simple computational alternative to general relativistic hydrodynamics , which 
can be used to probe the interaction of the torus with the black hole. Preliminary relativistic hydrodynamical 
simulations already reveal that this system gives rise to rich scenarios |4Q, i42] . 

• Stability of Kerr black holes under massive scalar fields. Although the stability question of Schwarzschild black 
holes has been settled the problem for spinning black holes is still quite open. In particular, it has been 
pointed out that for massive scalar fields, unstable modes are likely 49, 50]. A mathematical analysis of such a 
system is complicated even in the axisymmetric case (51| . Numerical simulations of massive scalar perturbations 
might help shed light on this problem. 

• Boson star orbiting a black hole [l^ l20l | . Such a simulation would shed light on both the orbital decay and the 
waveform radiated to null infinity. While this is perhaps the physically most interesting project, it is also the 
most computationally demanding because of the necessity to resolve the compact distribution of the scalar field. 

Some of the above projects can already be studied with the code as implemented in this work. However, for more 
"realistic" scenarios one would like to extend this treatment to include massive scalar fields and spinning black holes. 
Work in these directions is in progress. 
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TABLE I: Convergence in amplitude of the 3-dimensional EKG code 
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FIG. 1: The scalar field norm Q/X^ is plotted vs retarded time u for random initial scalar data with amplitude A = 10 * in a 
localized region between ra = 3 and = 5 on a grid of size Hx = 4:1, nq = Up = 25. The dashed line plots the rescaled norm 
Q/A^ for the same initial random data but with A = 2.5 x 10 ® on a grid of size Ux = 81, nq — rip = 45. 
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FIG. 2: The dashed Une plots the waveform r(j){u) (multiphed by 10^^) radiated to X+ in the q = p = direction resulting 
from the backscattering of an Z = 2, m = pulse. The pulse parameters for the initial data are ra = 3; rb = 5, A = 10^**. 
The grid size is 25 x 25 x 101. The solid line is the analytic curve corresponding to the complex quasi-normal mode frequency 
ui = 0.4875 — i0.098 obtained from a perturbation calculation. 
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FIG. 3: The rescaled waveform r4'{u)/\ is plotted for the q = p = Q direction at 1+. The parameters of the initial data are: 
= 3; r6 = 5, £ = 1, m = 0. The grid size is 25 x 25 x 101. All the curves overlap for A = 10"^, 10"*, 10"^. 



16 




17 




grid size is 89 x 89 x 101. 



18 




FIG. 6: Surface plots of the spin-weight invariant J, J on the North patch of T for u — 1,2, 3, 4, in zig-zag order from top to 
bottom, resuhing from the scattering of an axisymmetric pulse with initial data parameters ra = 3, rt = 5, A = 10~^. The grid 
size is 89 X 89 X 101. 




FIG. 7: Global radial-angular view of r<j) on the North patch for ii = 0, 2, 4, 6, 8, 10, in zig-zag order from top to bottom. We 
set p = 0, —1.2 < q < 1.2 (from bottom to top), 0.5 < a:: < 1.0 (from left to right). The initial data parameters are Va ~ 3, 
ri, = 5, A = 10"^ The grid size is 25 x 25 x 41. 
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FIG. 9: N+{u) (multiplied by 10^) is plotted for the g = p = 0.5 direction. The curves for A = 10"^, 10"^, 10"^ (broken line) 
overlap. For A = 10~^^ (continuous line), numerical error breaks the expected quadratic scaling. 



22 



1.5 




_2 I I I I I I I I I I I 

5 10 15 20 25 30 35 40 45 50 

u 

FIG. 10: The rescaled coefficient F2o{u)/\ (multiplied by 10^) is plotted for A = 10"^^, 10"^, 10"^, 10"^ all the curves overlap. 
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scattering of an asymmetric pulse. The initial data parameters are Va ~ 3.5, rj, = 12, A = 10 , fi = 0.03, qs ~ 0.2, ps — 0.3. 
The grid size is 85 x 85 x 101. 
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FIG. 14: Surface plots of the gravitational energy flux \N\'^ = |Af+|^ + \Nx\'^ on the South patch of for u — 1,2,3,4, in 
zig-zag order from top to bottom. The initial data parameters are ra = 3.5, rt — 12, A — 10"*^, fi = 0.03, qs ~ 0.2, ps = 0.3. 
The grid size is 85 x 85 x 101. The maximum value of |A''| is of order lO"'^^. 



